// display the convolved momentum distribution
// the input file format is mom str tot
// contributions from all configurations are displayed
// by Asia Sun, Jun. 23, 2022
// Last modified: Mar. 9, 2024

static int co = 0; // graph color
void primp(TGraph *g); // primp the graph
const string dir = "";
// N17-p c.s. on C12 target @240MeV/u
static const double csExp = 13.1, dcsExp = 1.5; // -p c.s.
static const char *cfg[] = {"0p10p1", "0p10p2", "0p32p1", "0p32p2"};
static const char *cfgltx[] = {"0p_{1/2}\\oplus0^{+}_{1}", "0p_{1/2}\\oplus0^{+}_{2}",
  "0p_{3/2}\\oplus2^{+}_{1}", "0p_{3/2}\\oplus2^{+}_{2}"};
static const double cmc = 17./16;
static const double ss[] = {0.76*cmc, 0.007*cmc, 0.08*cmc, 0.18*cmc};
static const int lsty[] = {2, 10, 3, 5}; // line style for each cfg

void momConv1(){
  string file = dir+"mom.txt";
  TGraph *gtot = new TGraph(file.c_str(), "%lg %*lg %lg"); primp(gtot);
  TLegend *lg = new TLegend(0.1290727,0.4461806,0.4849624,0.7517361,NULL,"brNDC");
  // TLegend *lg = new TLegend(0.6, 0.65, 0.89, 0.85);
  lg->SetBorderSize(0);
  lg->AddEntry(gtot, "Total");
  // add the configurations //
  const int nc = sizeof(cfg) / sizeof(const char *);
  vector<TGraph *> vg;
  char cfgfile[128]{}, gname[64]{};
  for(int i = 0; i < nc; i++){
    sprintf(cfgfile, "%s_totConv.txt", cfg[i]);
    vg.push_back(new TGraph(cfgfile, "%lg %lg"));
    sprintf(gname, "g%s", cfg[i]); vg[i]->SetName(gname);
    primp(vg[i]); vg[i]->SetLineStyle(lsty[i]);
    lg->AddEntry(vg[i], cfgltx[i]);
  } // end for over i

  TFile *f = new TFile("ppl.root");
  TH1F *hC16 = (TH1F *)f->Get("hmomdisNET"); hC16->SetLineColor(4);
  hC16->Scale(csExp/hC16->Integral("width")); // scale to the total c.s.
  lg->AddEntry(hC16, "^{17}N-p Mom. distri.");

  // scaling //
  // energy shift
  const double dm = gtot->GetMean() - hC16->GetMean() - 10;
  cout << "gtot->GetMean(): " << gtot->GetMean() << endl;
  cout << "hC16->GetMean(): " << hC16->GetMean() << endl;
  cout << "dm: " << dm << endl;
  // amplitude scaling
  const double cce = hC16->Integral("width")/gtot->Integral();
  // const double cce = 1.;
  cout << "gtot->Integral(): " << gtot->Integral() << endl;
  cout << "hC16->Integral(): " << hC16->Integral() << endl;
  cout << "cce: " << cce << endl;
  const int n = gtot->GetN(); double x, y;
  for(int i = 0; i < n; i++){
    gtot->GetPoint(i, x, y); gtot->SetPoint(i, x-dm, y*cce);
    for(int j = 0; j < nc; j++){
      vg[j]->GetPoint(i, x, y); vg[j]->SetPoint(i, x-dm, y*cce*ss[j]);
    } // end for over cfgs
  } // end for over i

	// const double dm1 = hC16p->GetMean() - hC16->GetMean();
	// cout << "dm1: " << dm1 << endl;
	// const double nx = hC16p->GetNbinsX();
	// for(int i = 1; i <= nx; i++) hC16p->SetBinContent(hC16p->-dm1, hC16p->GetBinContent(i));


  TCanvas *c = new TCanvas("c1", "c1", 800, 600);
  gtot->SetLineWidth(2); gtot->SetLineStyle(1);
  gtot->SetTitle("");
  gStyle->SetOptStat(0);
  TAxis *ax = hC16->GetXaxis(); ax->SetTitle("p_{c} [MeV/c]"); ax->CenterTitle();
  TAxis *ay = hC16->GetYaxis(); ay->SetTitle("d\\sigma_{-p}/dp_{c} [mb/(MeV/c)]");
  ay->CenterTitle(); ay->SetTitleOffset(1.30);
  hC16->SetTitle("");
  hC16->Draw();

  gtot->Draw("same");
  gStyle->SetOptStat(0);

  for(auto &g : vg) g->Draw("same");

  TF1 *gres = new TF1("gres", "gaus", 11600., 11900.);
  const double pc = hC16->GetMean(), maxc = hC16->GetMaximum();
  cout << "pc: " << pc << " maxc: " << maxc << endl;
  gres->SetParameter(0, maxc);
  gres->SetParameter(1, pc + 632.);
  gres->SetParameter(2, 0.003943*pc);
  lg->AddEntry(gres, "response function");

  lg->Draw();
  gres->Draw("same");
  c->Print("mom1.pdf");
  c->Print("mom1.root");
} // end main function mom

// primp the graph
void primp(TGraph *g){
  TAxis *ax = g->GetXaxis(); ax->SetTitle("p_{c} [MeV/c]"); ax->CenterTitle();
  TAxis *ay = g->GetYaxis(); ay->SetTitle("d\\sigma_{-p}/dp_{c} [mb/(MeV/c)]");
  ay->CenterTitle(); ay->SetTitleOffset(1.30);
  co++;
  if(5 == co || 10 == co) co++;
  g->SetLineColor(co); g->SetMarkerColor(co);
} // end function primp
